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Abstract We study by means of molecular dynamics 
simulations of periodic shear cells, the influence of par- 
ticle shape on the global mechanical behavior of dense 
granular media. Results at macro-mechanical level show 
that for large shear deformation samples with elongated 
particles, independent of their initial orientation, reach 
the same stationary value for both shear force and void 
ratio. At the micro-mechanical level the stress, the fabric 
and the inertia tensors of the particles are used to study 
the evolution of the media. In the case of isotropic parti- 
cles the direction of the principal axis of the fabric tensor 
is aligned with the one of the principal stress, while for 
elongated particles the fabric orientation is strongly de- 
pendent on the orientation of the particles. The shear 
band width is shown to depend on the particle shape 
due to the tendency of elongated particles to preferen- 
tial orientations and less rotation. 

Keywords Granular media • Anisotropy • Particle 
shape • Contact network ■ Shear cell 



1 Introduction 

, Granular materials present a complex mechanical response 
' when subjected to an external load. This global response 
l is strongly dependent on the discrete character of the 

medium. The shape, angularity and size distribution of 
^ the grains, the evolution of the granular skeleton (spatial 
, arrangement of particles, void ratio, fabric, force chains), 

and some phenomena occurring at the grain scale (like 

rolling or sliding) are determinant factors for the overall 

macroscopic response. 
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The study of anisotropy is of fundamental importance 
in order to understand the properties of granular ma- 
terials. In soil mechanics, Casagrande and Carrillo ij 
distinguished between inherent and induced anisotropy, 
the first as a result of the sedimentation of particles and 
the second as a product of inelastic deformation. Oda 
1 studied by means of optical microscopic photographs 
of thin sections of sand, the preferred orientations and 
spatial relation of the constituting grains, and found that 
independent of their shape particles tend to rest in a sta- 
ble configuration relative to the forces acting on them. 
Furthermore, he also found that the initial fabric is not 
only determined by the shape of the particles but also 
by the deposition or compaction method. 

In relation to inherent anisotropy Oda et al. Q and 
Oda and Nakayama have listed three sources in such 
materials: 

1. Anisotropic distribution of contacts (called structural 
anisotropy). 

2. Shape and preferred orientation of void spaces. 

3. Shape of the particles and preferred orientation of 
non-spherical ones. 

The complete alteration of inherent anisotropy due to 
types 1 and 2 during early stages of inelastic deformation 
in biaxial compression tests, on two-dimensional assem- 
blies of rods, was observed by Oda et al. They found, 
however, that the one due to type 3 was still present at 
large deformations. It would be therefore expected that 
in the so-called critical state of soil mechanics, which is 
associated with large shear deformation, the persistence 
of inherent anisotropy is mainly due to particle orienta- 
tion 5]. 

The role of anisotropy in granular materials has been 
widely investigated both experimentally and with nu- 
merical simulations. In relation to particle shape, a big 
effort has been recently done to properly characterize 
particle geometry and to study its effect on the me- 
chanical response of granular media. Concerning particle 
shape characterization Bowman et al. proposed a ex- 
perimental technique using fourier descriptors and image 
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Fig. 1 Schematic representation of a particle contact, the 
overlapping area a is indicated by the shaded zone. 



Up to now, however, there is no clear information at 
the micro-mechanical level about the influence of particle 
shape and the orientation of non-spherical particles on 
the evolution of granular materials and the correspond- 
ing anisotropic network of contacts towards the station- 
ary state reached at the global level under shear con- 
ditions. These subjects are very important in geotech- 
nical engineering and physics in order to get a better 
understanding of the mechanical response of granular 
materials. In this work, we use the molecular dynamics 
technique to simulate the mechanical behavior of two- 
dimensional particle assemblies contained in a periodic 
shear cell. The grains are represented by randomly gener- 
ated convex polygons. We focus on the influence of parti- 
cle shape on the overall plastic response. The dependency 
of the mechanical behavior on the evolution of inherent 
anisotropy (specially contact and non-spherical particle 
orientations) is studied. Results are analyzed from the 
macro and micro-mechanical point of view. 



analysis to assess particle morphology and texture, and 
in discrete element simulations (DEM) Matsushima and 
Saomoto 't'I presented a method to construct irregularly- 
shaped grains from shape data of real granular material. 
Regarding the influence of particle shape on the mechan- 
ical behavior of the granular media, Bowman and Soga 
1^ found that the stress-strain and creep response of fine 
silica sand is influenced by particle elongation or aspect 
ratio. Using two dimensional DEM simulations Shodja 
and Nezami and Nouguier-Lehon and Frossard 
studied the effect of particle elongation on the rolling and 
sliding mechanisms in granular media. Nouguier-Lehon 
et al. ^3 studied in biaxial test simulations the influ- 
ence of particle shape on the so-called critical state, and 
Peha et al. ^3 the effect on the asymptotical states of 
granular media. Three dimensional studies using ellip- 
soidal convex-shaped particles have also been performed 
highlighting the role of particle shape [Tsl : . 

The temporal and spatial evolution of anisotropic 
granular materials during cornpaction has also been in- 
vestigated. Villarruel et al. |l5j | found experimentally 
clear evidence that particle anisotropy can drive order- 
ing (configuration of particles evolves to a nematic state) . 
Lumay and Vandewalle jl6l | in experiments and simula- 
tions observed large variations of the asymptotic packing 
volume fraction as a function of the aspect ratio of the 
particles. Ribiere 0| using rice of different shapes found 
that grain anisotropy slows down the packing fraction 
evolution, and that the convection in the granular media 
and therefore the compaction mechanism depend also on 
it. 

The anisotropic network of contact forces in a gran- 
ular packing subjected to shear conditions has been re- 
cently investigated by Majmudar and Behringer [T^ in 
a two dimensional system of photoelastic discs and by 
Latzel et al. Ilfll in discrete element simulations. 



2 Model, sample construction and numerical 
experiment 

2.1 Molecular dynamics simulations 

We study dense polygonal packings, in which particles in- 
teract through visco-clastic contacts. The polygons can 
neither break nor deform, but they can overlap when 
they are pressed against each other. This overlap repre- 
sents the local deformation of the grains. This approach 
has been frequently used to model many different pro- 
cesses, such as strain localization and earthquakes ^20|. 
fragmentation , and damage . 

The repulsive, elastic normal contact force is propor- 
tional to the overlap area a between particles. In Figure 
n the configuration of a contact between two particles 
is presented. Pi and P2 represent the intersection points 
between the polygons; the segment that connects those 
points gives the contact line S = PiP2- This vector de- 
fines a coordinate system {n,t) at the contact, where 
t = S/|S| and n normal to it give the direction of the 
normal /„ and tangential ft components of the contact 
force. As point of application of the contact forces we 
choose the center of mass of the overlapping area [2^ . 
indicated in the figure with the dot inside the shaded 
zone. 

The elastic force at the contact point is calculated: 

r = -KSh - kS (1) 

where fc„ and kt are the normal and tangential contact 
stiffnesses, respectively. The deformation length 5 is cal- 
culated in terms of the overlapping area a and the length 
of the contact line |S|, 6 = a/|S|. The tangential friction 
force is introduced by an elastic spring whose length is 
equal to in the case of new contacts, while for old con- 
tacts the tangential spring- length is updated as follows, 



V^AtMD, 



(2) 
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where ^' is the previous length of the spring, AImd is 
the time step of the molecular dynamic simulation, and 
V( the tangential component of the relative velocity v° 
at the contact: 

^Vi- Vj -w, xl,+ Wj X Ij. (3) 

Vi is here the linear velocity and the angular ve- 
locity of the particles in contact. The tangential elastic 
displacement ^ at the contact may increase during the 
time that the condition \flr\ < /i/^ is satisfied. The slid- 
ing condition is enforced keeping constant the tangen- 
tial force when the Coulomb limit condition |/f | = /i/^ 
is reached. At this moment the permanent deformation 
begins. 

We also introduce a viscous force at each contact 
point, which is necessary to maintain the numerical sta- 
bility of the method and to take into account dissipation 
during particle's contact. 

^„ = -mK.v^r^^-^.,.v^^^) (4) 

where m = {l/rrii + l/mj)^^ is the effective mass (kg) of 
the two particles in contact, and and i^t are damping 
coefficients. Directly related to the restitution coefficients 
of the collisions . A background damping force of the 
form = —rnvbi'^i + w^) is also introduced, where i^b 
is the background damping coefficient. This background 
damping force is introduced in order to model the fric- 
tion between the particles and the bottom (or top) of 
the shear cell used on the two-dimensional experiments 
performed by Veje et. al and Howell et al. p5l |. 



2.2 Generation of samples 

In order to study the influence of particle shape on the 
mechanical behavior of granular media, two different ty- 
pes of convex polygons were used to perform the numer- 
ical simulations. Figure [51 shows the visual differences 
between the two types of particles. Polygons depicted in 
Figure Ela) with an almost isotropic shape, from now 
on, will be called 'isotropic' particles, and polygons in 
El^b) will be referred as 'elongated' particles. The shape 
of the particles will be characterized by the aspect ratio 
A, which is obtained from the ratio between the length 
of the longest and shortest axis of the particles. 

The random generation of the polygons is accom- 
plished by means of a Voronoi tessellation . A regular 
square lattice is used to generate the isotropic polygons. 
This is done by setting a random point in each cell of 
the lattice. The polygons are constructed by assigning 
to each point the part of the plane that is nearer to it 
than to any other point. After generation, particles are 
moved apart by multiplying their coordinates by a con- 
stant in order to obtain a very loose state. Particles are 
then compressed isotropically by four rigid walls until 
the desired conflning pressure is reached. The elongated 
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Fig. 2 Two types of particles used in the numerical simula- 
tions: (a) isotropic A = 1.0 and (b) elongated polygons A = 
2.3. Points are the center of mass of the polygons. 

particles are generated by stretching or contracting, be- 
fore the compression, the particles in horizontal or/and 
vertical direction by a given factor. The average elon- 
gation of the grains is given by the ratio between these 
factors A (aspect ratio) . In our numerical simulations we 
consider two aspect ratios: A = 1.0 for isotropic and A 
= 2.3 for elongated polygons (which are clearly anisotro- 
pic). Although aspect ratios larger than 2.3 can be used, 
we choose this value in order to avoid particles with very 
sharp angles that could have an unrealistic overlap. Nev- 
ertheless, this should be verified by means of a system- 
atic study of different and/or larger aspect ratios of the 
particles. Finally, by setting the interparticle friction to 
zero during the compression we obtain the dense samples 
used in our simulations. 



2.3 Numerical experiment 

Wc want to study the behavior of the granular media 
in a periodic shear cell. The cell contains 1500 particles, 
being 50 particle diameters wide and 30 diameters high. 
Periodic boundary conditions are imposed in horizon- 
tal direction. The top and bottom have fixed boundary 
conditions. A constant confining stress (or normal force) 
is imposed between the bottom and the top horizontal 
walls. The top and bottom layers of particles are sheared 
in opposite direction with a fixed horizontal velocity 
relative to each other. The particles in one layer are not 
allowed to rotate or move against each other. The top 
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Fig. 3 Sketch of the shear cell. A normal force is applied between top and bottom wall. A horizontal velocity is fixed 
shearing the sample. Green (grey) particles correspond to the image used to implement the periodic boundary conditions. 



boundary is free to move in vertical direction in order 
to permit a volumetric change of the sample, while the 
bottom is kept fixed. A configuration of the shear cell 
used in our simulations is depicted in Figure |3| 

In all simulations the mechanical parameters of the 
particles are: interparticle friction coefficient /i = 0.5, 
normal stiffness at the contact fc„ — 1.6 • 10* N/m, and 
normal damping coefficient Vn = 4000 s^^. The stiffness 
ratio C, = kt/kn as well as the viscosity ratio vt/i^n were 
taken to be 1/3. A background damping coefficient v\, 
= 12 s""'^ was used. Further simulations, results not pre- 
sented here, with no and different background damping 
coefficients z/f, in order to evaluate its influence on the 
mechanical behavior of the medium were performed. We 
found that the damping Vf, in the range used here has 
no effect on the evolution of the internal variables of the 
medium, and additionally that the Vb value used in this 
work has no influence on the global mechanical response 
either. 

The horizontal and vertical directions are indicated 
as X and y, respectively. In order to study the evolution of 
the packing we use the strain variable 7, which is defined 
as follows: 

7 = D,/ho, (5) 

where is the horizontal displacement of the boundary 
particles and ho is the initial height of the sample. The 
void ratio e of the sample is related to the volumetric 
deformation, and is calculated: 

e = VT/Vs-l, (6) 

where Vr is the total volume of the sample and Vs the 
volume occupied by all the particles. 



3 Shear test simulations 

3.1 Global mechanical behavior (effect of initial 
configuration) . 

3.1.1 Statistically different samples 

Samples corresponding to different seeds for the random 
number generation of the Voronoi tesellation are used to 
evaluate the global mechanical response of the granular 
packing. This is done in order to assess whether different 
initial configurations of particles reach the same steady 
state. Our results correspond to a shear velocity of 10 
cm/s, and to elongated polygons initially oriented per- 
pendicular to shear direction. In Figure^jthe evolution of 
the resultant shear force and the void ratio is presented 
for the different configurations. In Figure 01a), the shear 
force Fs is normalized by the normal force F„ applied to 
the system. Initially, the ratio Fs/Fn has a strong incre- 
ment related to the breaking of the interlocking of the 
particles. After this stage, a saturation towards a nearly 
constant value of the Fs/Fn ratio necessary to shear the 
granular media is observed, this behavior is identical for 
all the samples. For small values of strain, the evolution 
of the void ratio (Figure ^p) also presents a high initial 
increase saturating later at a constant value. This satu- 
ration occurs slower than for Fs/Fn. Samples with elon- 
gated particles saturate at a higher value of Fg/Fn and 
also higher void ratio, as a consequence of the stronger 
interlocking due to the particle shape. One can there- 
fore conclude that elongated grains are more sensitive to 
volumetric changes and develop a higher shear strength. 
This result had been in fact previously observed [3| ■ 
We consider this saturation of the Fs/Fn value and void 
ratio e as the steady state of the sheared material. 
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Fig. 4 Evolution of (a) sliear force and (b) void ratio 
for different samples with the same mechanical parameters. 
Isotropic (A = 1.0) and elongated particles (A — 2.3) are rep- 
resented by solid and dashed lines, respectively. 
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Fig. 6 Evolution of (a) shear force and (b) void ratio for 
samples with different initial particle orientation, A = 2.3. 
Samples labelled H and V have particles oriented in horizon- 
tal and vertical direction, respectively. R corresponds to an 
initial random rotation of the particles. 
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Fig. 5 Evolution of the coordination number for different 
samples with the same mechanical parameters. Isotropic par- 
ticles A = 1.0 (solid lines) and elongated ones A = 2.3 (dashed 
lines) . 



The evolution of the coordination number for isotropic 
and elongated particles is depicted in FigureEl Note that, 
despite of reaching a higher void ratio, samples with elon- 
gated particles saturate at larger value of coordination 
number compared to isotropic particles. This can be un- 



derstood in terms of a geometrical effect, consequence 
of the flat shape and/or larger relative plane surface in 
the A = 2.3 case, which allows for a higher number of 
contacts per particle. 

3.1.2 Different initial particle orientations 

We study in this Section the influence of the anisotropy 
on the macroscopic behavior of granular media due to the 
initial orientation of elongated particles. Three different 
initial configurations are obtained for the samples used 
in this analysis: 

1. On average the grains are oriented parallel to the 
shear direction (this will be called "horizontal" sam- 
ple - H). 

2. On average the grains are oriented perpendicular to 
the shear direction (we will call this the "vertical" 
sample - V). 

3. Grains (H or V) are randomly rotated before isotropic 
compression (which we call the "random" samples - 
HR or VR). 

Configurations number 1 and 2 correspond to samples 
with different initial orientation of the particles. Config- 
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uration number 3 is equal to configurations 1 and 2, but 
with an additional induced random rotation to the par- 
ticles before the compression (between and 27r rad). 
In all three cases the samples wind up having a slight 
deviation from the originally induced anisotropy due to 
particle interactions during the compression. 

In Figure El the evolution of Fg/Fn and the void ra- 
tio e for samples with particles initially oriented in hori- 
zontal and vertical direction, and an additional random 
rotation is presented. Results correspond to a shear ve- 
locity of 40 cm/s. We notice that Fg/Fn and the void 
ratio e evolve toward the same saturation value when 
they reach the steady state independently of the initial 
anisotropy due to contact and particle orientations. This 
independence of the initial anisotropy will be explained 
by studying the evolution of the internal variables in Sec- 
tion E3 



3.2 Evolution of internal variables 

In order to get a better understanding of what is occur- 
ring at the micro-mechanical level and to find some ex- 
planations for the macro-mechanical behavior observed 
in the previous section, we study the evolution of the lo- 
cal stress, the fabric and the inertia tensors of the sam- 
ples. Three samples are used: one with isotropic poly- 
gons, and two with elongated polygons (corresponding 
to an initial horizontal and vertical orientation of the 
particles, as explained in the previous section). Before 
showing these results, some key concepts that will be 
required for the description of the media are introduced. 



3.2.1 Definitions 

The anisotropy of the contact network within the gran- 
ular sample can be characterized by the fabric tensor of 
second order F. The fabric tensor takes into account the 
distribution of the orientations of the contacts between 
particles, that is to say the geometrical structure of the 
medium Q. For a single particle p its components Fij 
are obtained from 



pp. 



where the dyadic product of the so-called branch vector 
Z'^, connecting the center of mass of the particle to the 
contact point c, is summed over all the contacts Cp of 
particle p. Note that, so defined, the trace of the fabric 
tensor Ff^ is then the number of contacts Cp of the par- 
ticle p. It is also possible to define a normalized fabric 
tensor F[j/Cp, whose trace is unity. Finally, the mean 
fabric tensor for an assembly can then be defined as fol- 
lows: 



F,., 



1 



(8) 



where the particle fabric tensor F[j is summed over the 
total number of particles Np within a representative vol- 
ume element (RVE). The trace of this tensor is the lo- 
cal mean coordination number Cm, and therefore the 
normalized mean fabric tensor can also be defined as 

Fij I CjYi . 

The inertia tensor can be calculated for each particle 
as follows: 



j p{5ij E^fc ~ XiXj)dA 



(9) 



where p is the density of the particles (kg/m^), 5ij is the 
Kronecker delta symbol, k runs in our two-dimensional 
case from 1 to 2, dA is the differential area element, and 
X is the shortest distance from the rotation axis to dA. 
The mean inertia tensor is calculated in terms of ©: 



E^ 

p=i 



(10) 



where the particle inertia tensor i^^ is summed over the 
total number of particles Np in the RVE. 

We also calculate the stress tensor, which is defined in 
terms of the contact force /J between the grains (acting 
at the contact point c), and the branch vector belong- 
ing to the contact point. The stress tensor is defined as 
follows W\: 



1 



(11) 



where V is the volume of the RVE, and the summatory 
extends over all the contacts Nc in the RVE. 

Finally, we will compute the principal directions (ma- 
jor M and minor m) of the mean fabric F, inertia I and 
stress a tensors. These principal directions are measured 
with the horizontal axis x. We denote Op the major prin- 
cipal direction of the fabric tensor, Oj the major princi- 
pal direction of the inertia tensor and 9a the one of the 
stress tensor. We are also interested in the individual ori- 
entation of non-isotropic particles. We designate 9^ the 
angle made by the major principal direction of the iner- 
(7) tia tensor i^^ of the particle p with the horizontal axis 
X. Therefore, 9^ gives the preferred orientation for each 
particle of the assembly 

3.2.2 Evolution of the internal variables 

In Figure Efa-c), we show the orientational distribution 
of the branch vectors and the principal directions of the 
mean fabric tensor {Fm and F„i) for the initial configura- 
tion of the samples. Observe that, in the case of isotropic 
polygons (Fig. El^a)), the distribution presents no pre- 
ferred direction within the statistical fluctuations. For 
elongated polygons (Fig. [7fb-c)), however, one can ob- 
serve that the major principal component of the fabric 
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Fig. 7 Polar distribution of branch vectors in the initial configuration (a,b,c) and the steady state (d,c,f), for isotropic 
particles A = 1.0 (a,d), and elongated particles (A — 2.3) initially oriented in horizontal direction (b,e) and in vertical 
direction (c,f). The principaJ directions of the mean fabric tensor {Fm and Fm), and the reference axes x and y are plotted 
with solid and dashed lines, respectively. The radius of the dashed circle corresponds to the maximum value of the distribution. 
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Fig. 8 Polar distribution of particle orientations 6^, initial configuration (a,c) and in the steady state (b,d) for elongated 
particles (A = 2.3) initially oriented in horizontal direction (a,b) and in vertical direction (c,d). The principal directions of 
the global inertia tensor (7m and Im), and the reference axes x and y are plotted with solid and dashed lines, respectively. 
The radius of the dashed circle corresponds to the maximum value of the distribution 
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Fig. 9 Evolution of the deviatoric component of the fabric 
(a) and the inertia tensor (b) for isotropic particles A = 1.0, 
and elongated particles (A = 2.3) initially oriented in hori- 
zontal direction (H) and in vertical direction (V). 
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Fig. 10 Evolution of the quotient of the principal compo- 
nents of the fabric (a) and the inertia (b) tensors for isotropic 
particles A = 1.0, and elongated particles (A = 2.3) initially 
oriented in horizontal direction (H) and in vertical direction 
(V). 



tensor F is oriented towards the direction in which the 
polygons were initially stretched. 

The anisotropic distribution of the contact orienta- 
tions is more pronounced in case (b) than in case (c); 
this is most probably due to the shape of the shear cell 
(which indeed is wider than higher) and the compression 
process using rigid walls. The orientational distribution 
of the contacts in the steady state is depicted in Fig- 
ure EJd-f). We notice that the distribution of contact 
orientations for elongated particles (Fig. EJe-f)) is very 
similar independent of their initial orientation, while for 
isotropic particles (Fig. I3d)) it is clearly different. The 
major principal direction of the fabric tensor follows this 
same trend. 

In FigurelHl the polar distribution of 6"' for elongated 
particles, and the principal directions of the mean in- 
ertia tensor (/m and Im) in the beginning and in the 
stationary state are presented. We observe that, similar 
to the case of contact orientations, they evolve towards 
the same global orientation independently of the initial 
particle directions. 

In order to study the evolution of the fabric and 
inertia tensors we monitor their deviatoric component 
Fyy — Fxx, and the quotient F^x / Fyy during the simu- 
lation. One can observe in Figure El where the evolution 



of the deviatoric component of Fij and lij is shown, that 
the deviatoric reaches a stationary value for both types of 
particles, and that the induced anisotropy is much higher 
for elongated particles than for isotropic ones. The same 
result is observed for the quotient of the principal com- 
ponents of the tensors (FigureEHl- This stationary value 
of the deviatoric component and the quotient is directly 
related to the steady state at the macro-mechanical level, 
and seems to be a micro-mechanical requirement for the 
global steady state. This assumption is supported by 
simulations of biaxial tests reported by Nouguier-Lehon 
et al. 1^ , where samples with elongated particles do not 
reach neither a stationary value for the components of 
the fabric and the orientation tensors nor the so-called 
critical state at the macro-mechanical level, but samples 
that reach the stationary state for the components of the 
tensors do so at the global level. 

Furthermore, for samples with elongated polygons 
the deviatoric part of the Fij and /y tensors, and the 
ratio of their principal components reach approximately 
the same stationary value independent of the initial par- 
ticle orientations. This means that the initial inherent 
anisotropy (fabric and particle orientation) is completely 
erased and reoriented in direction of the induced shear 
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Fig. 11 Evolution of the major principal direction of the fab- 
ric (a), the inertia (b), and the stress (c) tensors for isotropic 
particles A = 1.0, and elongated particles (A — 2.3) initially 
oriented in horizontal direction (fl) and in vertical direction 
(V). 



during the experiment. The evolution of the major prin- 
cipal direction of the fabric 9f, the inertia Oj and the 
stress tensors 9a are shown in Figure ^3 In the case of 
the inertia tensor, the major principal direction 6j is re- 
oriented for all samples towards an angle close to 160°. 
For the stress tensor, 9a- (~ 45°) is the same for both 
particles independent of the particle shape. This orien- 
tation of the stress comes from the direction of the force 
chains carrying the largest stresses (Fig.^Jc-d)). On the 
other hand, the major principal direction of the fabric 
tensor ^i;^ is completely different for isotropic and elon- 



gated polygons. For isotropic particles it is similar to the 
major principal direction of stress 9a-, but for elongated 
polygons the fabric orientation is close to the one of the 
inertia tensor. 

To clarify this result further, the contact forces larger 
than the mean value and the principal axes of the fab- 
ric tensor of the corresponding particles are plotted in 
Figure El This is done for both types of polygons, and 
for the initial configuration and a snapshot in the steady 
state. In Figure ITW a'). where the initial configuration of 
isotropic particles is shown, one can observe that the ma- 
jor principal axis of the fabric tensor of each particle F^j 
is oriented independently of the orientation of the force 
chains. In the system with elongated particles, however, 
is slightly oriented in the largest dimension of the 
particles (major principal axis of the tensor of inertia of 
each particle i^/). 

In the stationary state, we notice that in the system 
with isotropic particles F^j approximately follows the 
direction of the force chains that carry the larger forces. 
In the one with elongated particles, on the contrary, F^j 
is oriented in direction of the largest dimension of the 
particles i^j . The orientation of the elongated particles 
within the main force chains is associated to the stability 
of the packing. That is to say; forces are transmitted 
through contacts closely parallel to the minor principal 
axis of the inertia tensor of each particle zJJ^, in which 
flat contact surfaces give a more stable configuration to 
the system. We conclude then that the orientation of the 
contacts in the steady state, in the case of non-spherical 
particles is governed by the particle orientation, and for 
isotropic particles by the direction of the major principal 
stress. This is also observed in Figure ITTl where in the 
steady state for isotropic particles 9f is almost the same 
as 9a, and for elongated particles 9f is nearly 9j. 



3.3 Shear localization and particle rotation 

In order to study strain localization and particle rota- 
tion, the shear cell is divided into horizontal layers, i.e 
parallel to shear direction. For a clearer presentation of 
the results, we normalize the vertical dimensions with 
the height of the system h. The origin corresponds to 
the bottom and 1 to the top of the sample. We use in 
our analysis the rotation that particles accumulate dur- 
ing every unit increment (Ajunit) of the strain variable 
7 in the steady state (in our experiment we take j-miuai 
= 10 and Ifinai — 35, i.e. in total 25 Ajunit)- Then, we 
average this accumulated particle rotation for each layer 
of the system, and for all the considered strain incre- 
ments. In Figure 1131 the average accumulated rotation 
in the steady state for each layer and a shear velocity of 
40 cm/s for isotropic and elongated particles initially ori- 
ented in vertical direction is shown. We observe a clear 
localization of rotations, having a peak close to the center 
and decreasing toward the boundaries. This distribution 
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Fig. 12 Force chains (light hnes, thickness proportional to magnitude) and principal axes of the fabric tensor (black lines), 
for initial configuration (a,b) and the steady state (c,d), for isotropic A = 1.0 (a,c) and elongated particles A = 2.3 (b,d). 



resembles the movement of two rigid bodies against each 
other on a shear band. 

We are also interested in the width of the shear zone, 
which we define here as the width of the zone with parti- 
cle rotation larger than 80 % of the maximum rotation. 
One can notice that there are two important differences 
between the two types of particles considered: 

— The average accmnulated rotation is lower for elon- 
gated particles, in this particular case it is 65 % of 
the isotropic particles rotation. 

— The width of the localization zone is smaller for elon- 
gated particles (0.45 times the system height h for 
elongated and 0.55 times h for isotropic particles). 

These differences in accumulated rotation and rel- 
ative width between isotropic and elongated particles 
can be explained by the frustration of movement and 
rotation that elongated particles experience due to the 
stronger interlocking among them. In this way, the lo- 
calization zone (rotation zone) for elongated polygons 
becomes thinner than for isotropic ones. In Figure 1141 
the mean accumulated rotation since the beginning of 
the simulation for isotropic and elongated particles is 
depicted. Notice that the mean rotation is almost twice 
for isotropic than for elongated particles at the end of 
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Fig. 13 Average accumulated rotation of the particles dur- 
ing the steady state, isotropic A = 1.0 (full dots) and elon- 
gated particles A — 2.3 V (open squares) within horizontal 
layers as a function of relative depth. The maximum and 80 
% of the maximum value of rotation are plotted as solid and 
dashed lines, respectively. 



the simulation. We also calculate the probability distri- 
bution function of particle rotation, which is shown in 
Figure [Tsl for both isotropic and elongated particles and 
7 = 35. For isotropic particles a more uniform distri- 
bution is observed, and the maximum value is close to 
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Fig. 14 Mean accumulated rotation since the beginning of 
the simulation for isotropic particles A = 1.0 (light line), and 
elongated particles (A = 2.3) initially oriented in horizontal 
direction (H, black line), and in vertical direction (V, black- 
dashed line). 




Fig. 15 Probability distribution function of accumulated ro- 
tation since the beginning of the simulation until 7 — 35, for 
isotropic particles (A — 1.0, light line), and elongated parti- 
cles (A = 2.3) initially oriented in horizontal direction (black 
line), and in vertical direction (black-dashed line). 



four complete rotations (a complete rotation 2t: rad) . For 
elongated particles the probability distribution function 
presents several peaks every half of rotation (tt rad) . This 
fact indicates the strong frustration of rotation that such 
particles undergo during shearing, and that the typical 
mode of accumulating rotation is then every half com- 
plete rotation. 



4 Discussion and final remarks 

In this paper, the influence of particle shape on the me- 
chanical behavior of sheared dense granular media has 
been investigated using two-dimensional discrete numer- 
ical simulations. Our results show an important influ- 
ence of particle anisotropy on both the macro and micro- 
mechanical behavior of the granular media. It is found 
that for samples with isotropic and elongated particles 



the shear force and volumetric strain saturate at constant 
values reaching a steady state. These values in the case 
of elongated particles are higher than for isotropic parti- 
cles due to the stronger interlocking between anisotropic 
particles. This is consequence of the geometrical effect 
of particle shape. Furthermore, samples with elongated 
particles reach the same saturation value in the steady 
state independently of the initial orientation of the par- 
ticles. This is related to the removal and reorientation of 
the initial inherent anisotropy (fabric and particle orien- 
tations) in direction of the induced shear. This has been 
confirmed by studying the evolution of the fabric and the 
inertia tensors. 

The deviatoric part and the quotient of the principal 
components of the fabric F and inertia I tensors, for both 
types of particles, reach a stationary value independent 
of their initial one. This is directly related to the steady 
state at the macro-mechanical level. The principal direc- 
tions of F and I present the following behavior: in the 
initial state of the samples, isotropic polygons present 
no preferred direction of contacts, however, in the case 
of elongated polygons the major principal component of 
F is oriented along the direction of the major principal 
component of I. In the steady state, Fm in the case of 
isotropic particles is reoriented in the direction of the ma- 
jor component of the stress tensor a, but for elongated 
particles Fm evolves following the induced orientation 
that particles undergo during shearing. The direction of 
the major component of a is the same for both particle 
shapes. Independently of the initial orientation, samples 
with elongated particles reach the same contact Op and 
particle Oj global orientation in the steady state. One 
can then conclude that a stationary value of the princi- 
pal components and principal directions of the fabric and 
inertia tensors is a micro-mechanical requirement for the 
existence of the global steady state of the medium. We 
also conclude that for the sheared dense granular mate- 
rials here presented, the contact orientation in the global 
stationary state is governed in the case of isotropic par- 
ticles by the direction of the major principal component 
of the stress tensor, and for elongated particles mainly 
by the major principal component of the inertia tensor 
(particle orientation) . 

At the particle level, these results are clearly under- 
stood when the inertia and fabric tensor of the parti- 
cles within the force chains carrying the larger forces are 
studied. We found that the orientation of elongated par- 
ticles is associated to the stability of the packing, i.e. 
forces are transmitted through contacts parallel to the 
shortest dimension of the particles i^. 

Concerning strain localization and particle rotation 
the following is observed: the width of the shear zone and 
the accumulated rotation is larger for isotropic particles 
than for elongated particles. This results is explained by 
the frustration of rotation that elongated particles expe- 
rience due to the stronger interlocking among them. This 
is clearly observed in the probability distribution func- 
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tion of particle rotation, where the typical mode of accu- 
mulating rotation for elongated particles is every tt rad. 

Results presented in this paper lead to a better com- 
prehension of the role of particle shape on the macro 
and micro-mechanical response of granular materials. It 
is, however, necessary to say that due to the nature of 
our two-dimensional analysis, a validation with a three- 
dimensional model would be very helpful. Further work 
will focus on a systematic study of the influence of differ- 
ent aspect ratios of the particles, beyond the two limit 
cases here studied, on the observed phenomena. Addi- 
tionally, the influence of two types of particles (isotropic 
and anisotropic) as constituents of the sample on the 
mechanical behavior of the granular media will be inves- 
tigated. 
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